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ABSTRACT 

Context. The continuum intensity at wavelengths around 1 mm provides an excellent way to probe the solar chromosphere and thus 
valuable input for the ongoing controversy on the thermal structure and the dynamics of this layer. 

Aims. The synthetic continuum intensity maps for near-millimetre wavelengths presented here demonstrate the potential of future 
observations of the small-scale structure and dynamics of internetwork regions on the Sun. 

Methods. The synthetic intensity / brightness temperature maps are calculated on basis of three-dimensional radiation 
(magneto-)hydrodynamic (MHD) simulations. The assumption of local thermodynamic equilibrium (LTE) is valid for the source 
function. The electron densities are also treated in LTE for most maps but also in non-LTE for a representative model snapshot. 
Quantities like intensity contrast, intensity contribution functions, spatial and temporal scales are analysed in dependence on wave- 
length and heliocentric angle. 

Results. While the millimetre continuum at 0.3 mm originates mainly from the upper photosphere, the longer wavelengths consid- 
ered here map the low and middle chromosphere. The effective formation height increases generally with wavelength and also from 
disk-centre towards the solar limb. The average intensity contribution functions are usually rather broad and in some cases they are 
even double-peaked as there are contributions from hot shock waves and cool post-shock regions in the model chromosphere. The 
resulting shock-induced thermal structure translates to filamentary brightenings and fainter regions in between. Taking into account 
the deviations from ionisation equilibrium for hydrogen gives a less strong variation of the electron density and with it of the optical 
depth. The result is a narrower formation height range although the intensity maps still are characterised by a highly complex pattern. 
The average brightness temperature increases with wavelength and towards the limb although the wavelength-dependence is reversed 
for the MHD model and the NLTE brightness temperature maps. The relative contrast depends on wavelength in the same way as 
the average intensity but decreases towards the limb. The dependence of the brightness temperature distribution on wavelength and 
disk-position can be explained with the differences in formation height and the variation of temperature fluctuations with height in the 
model atmospheres. The related spatial and temporal scales of the chromospheric pattern should be accessible by future instruments. 
Conclusions. Future high-resolution millimetre arrays, such as the Atacama Large Millimeter Array (ALMA), will be capable of 
directly mapping the thermal structure of the solar chromosphere. Simultaneous observations at different wavelengths could be ex- 
ploited for a tomography of the chromosphere, mapping its three-dimensional structure, and also for tracking shock waves. The new 
generation of millimetre arrays will be thus of great value for understanding the dynamics and structure of the solar atmosphere. 

Key words. Sun: chromosphere, radio radiation; Submillimeter; Hydrodynamics; Radiative transfer 



1. Introduction 

Many details concerning the small-scale structure of the solar 
chromosphere remain an open issue despite the large progress 
during the last decades on the observational but also on the 
modelling side. Recent works suggest that the solar chromo- 
sphere within internetwork regions is a highly inhomogeneous 
and dynamic phemenon demandin g for high spatial and tempo- 
ral re solution on bo th sides (e.g.. lAl et al.l 120021: i Krii ger et al J 
200 It l^ res 200| IWoger et alJ |2006; TritschleiUdJ |2007t 
Vecchio et al.il2007 . and references therein) On the numerical 
side, the increasing computational power but also sophisticated 
new methods now allow for the necessary resolution. On the 
other hand, advanced instruments enable new highly-resolved 
observations of hitherto not achieved quality (e.g., with SOT 
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onboard the Hinode satellite, see, e.g., iTsunetal 12006*). For a 
detailed comparison of observed images and three-dimensional 
simulations, synthetic intensity maps need to be calculated. 
Unfortunately, the chromospheric diagnostics so far available, 
such as the calcium resonance lines, require the treatment of de- 
viations from local thermodynamic equilibrium (LTE). A proper 
treatment demands for large computational resources and so- 
phisticated methods, rendering such calculations hardly feasi- 
ble for 3D models. And even in case of successful calculations, 
the complex translation of temperature into emergent intensity 
complicates the interpretation and hampers deriving the thermal 
structure from observations. 

An exception, in contrast to most other diagnostics, is the 
radio continuum at millimetre and sub-millimetre wavelengths. 
As its source function can be treated in LTE, it can be syn- 
thesised easily and hence offers a co nvenient way to c o mpare 
observations and numerical models. iLoukitcheva et al.l (|2004 
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I2OO6I) took advantage of this fact and compared a large collec- 
tion of observations in the millimetre and sub-millimetre range 
wi th synthetic brightnes s temperatures calculated from models 
bv iFontenla et al.l1l993L herea fter FAL) and [C arlsson & Steiiil 
(I1995L hereafter CS). Recently IToukitcheva et al., (i2006l) com- 
pared these models to observat ions done with th e Berkeley- 
lUinois-Maryland Array (BIMA, IWhite et alJ l2006V The major 
problem of observations at millimetre wavelengths is the gen- 
erally poor spatial resolution which renders granular scales so 
far unaccessible - even with the BIMA array with its 10 anten- 
nae. The situation will substantially improve with the next gen- 
eration of large interferometric arrays, e.g., the Atacama Large 
Millimeter Array (ALMA ) which will be fully operational in 
~ 2012 (e.g. iBeaslev et al. l l2006l) . This instrument will provide 
high spatial and temporal resolution, finally allowing to observe 
the small-scale structure of the solar chromosphere in detail. 

Here we use three-dimensional radiation hydrodynamics 
simulations to synthesise the continuu m intensity at millimetre 
wave lengths from 0.3 mm to 9 mm (see lWedemever-Bohm et aU 
12005 1 for a precursory study). Since the simulations do only in- 
clude weak magnetic fields or no field, the present analysis refers 
to internetwork regions only. 

In Sect. |2] some instruments, that are potentially interesting 
for solar observations, are introduced. The numerical simula- 
tions and the method of producing synthetic intensity images 
are described in Sects. [3] and H] respectively. The results are 
presented in Sect. |5] followed by discussion and conclusions in 
Sect. |6] and Sect.|2l 



olution Aa we assume full u-v coverage in the limit of a large 
number of baselines. We finally derive the relation 



(Aa), 



Vm VA^a(A^a-l) 



(2) 



The primary beam size d and thus the resolution, too, depends 
proportionally on wavelength A. Note that our estimate is not 
a hard limit. Reliable images should be possible at signifi- 
cantly higher resolution, using the known positivity of the im- 
age and other 'a priori' constraints. The technique of Multi- 
Frequency Synthesis (MFS) uses the frequency-dependence of 
the u-v coverage, re sulting in a larger number of data constraints 
(IConw ay et al.l[T99 0). It has thus the potential of significantly 
increasing the effective resolution of reliable images. However, 
detailed imaging simul ations would be needed to quantify the 
achievable resolution ( Conwavll2007h . 

The list of interferometers in Sects. ITl 112. 5l contains charac- 
teristic properties that are relevant for the discussion in Sect. 16.61 
Some properties like the FOV, however, are not easily expressed 
with a single number - in particular for the heterogeneous ar- 
rays CARMA, FASR, and RAINBOW since the different dish 
diameters of the individual antenna types result in different 
(wavelength-dependent) primary beam sizes. Please note that the 
description of the interferometers relies on information available 
in the literature and in the internet and is only thought to give a 
broad overview. Technical details might differ slightly in prac- 
tice. 



2. Instruments 

The instruments addressed in this section can (optionally) be 
used as interferometers and are potentially interesting for solar 
observations. Although the angular resolution of an interferome- 
ter can be very large depending on the maximum baseline (con- 
necting two individual antennae), there is an effective maximum 
resolution on which highly reliable images of objects such as 
the Sun, with complicated emission filling the whole primary 
beam of each antenna, can be made. Being a dynamic object the 
amount of received data constraining a snapshot image equals 
twice the number of baselines (real and imaginary parts of each 
visibility being counted separately). The number of unknowns 
for the image construction, however, is equal to the number of 
independent synthesised beam areas within the primary beam. 
Reliable imaging requires that there are more data constraints 
than there are unknowns. This condition defines the effective 
resolution of reliable images for the interferometer. Here we es- 
timate this resolution by assuming a homogeneous distribution 
of the synthesised beams ("image elements") over the primary 
beam area (the field of view, FOV). The number of those el- 
ements is equal to the number of baselines. An array with A^a 
antennae gives a maximum of 

iv. = M|zi) „, 

possible basehnes. Depending on technical details of the array, 
only a subset might be realised. As such details might change, 
we always refer to the maximum number in this work. Also note 
that a necessarily finite number of base lines limits the number 
of resolution elements. A finite number of baselines, however, 
can only provide an incomplete coverage of the M-v-plane (spa- 
tial Fourier space), making it difficult to determine the effective 
spatial resolution. For our estimate of the effective spatial res- 



2.1. ALMA 

Among a large range of issues important to modern astronomy, 
the Atacama Large Millimeter Array (ALMA) will also be used 
for observations of the Sun. While its single-dish predecessor 
APEX (Atacama Pathfinder Experiment) has already been 
installed and delivers first scientific results, the construction of 
the array started next to the APEX site on a plateau at 5000 m al- 
titude in the Chilean Andes. The science verification phase will 
most likely start in late 2009, followed by an early science stage 
from 2010, and finally full operation in 2012. All details given 
in this section refer toBastian ( 2002h. lBrown etalJ (|2004), more 
recently IBeaslev et al.l 12006 ) and the ALMA web pages (e.g., 
http : //www . alma . irif o, http : //www . alma . itrao . edu, 
http : //www ■ esc . q rg/pro jects/alma/). See also 
lEscoffieret all (120071) . 

The main array will consist of 50 antennae with a diame- 
ter of 12 m in an adjustable configuration ranging from a com- 
pact size of 150 m to a maximum base length of ~ 18.5 km. It 
is supplemented with the (semi-independent) Atacama Compact 
Array (AC A) with 12 antennae with 7 m diameter and 4 12 m- 
antennae. Each ALMA antenna will be equipped with receivers 
which cover up to ten frequency bands. Initially only the fre- 
quency bands in the range from 84 GHz to 373 GHz will op- 
Table 1. Frequency bands of ALMA that will be installed first 
and corresponding values for the expected primary beam size 
diameter d,i and the estimated spatial resolution (Aa), for images 
of maximum reliability (see text for details). A total number of 
50 antennae is assumed. 



band 


V [GHz] 


A [mm] 


dA ["] 


(Aa),, ["] 


3 


84 - 116 


3.57 - 2.58 


75 - 54 


1.51 - 1.10 


6 


211 - 275 


1.42 - 1.09 


30 - 23 


0.60 - 0.46 


7 


275 - 373 


1.09 - 0.80 


23 - 17 


0.46 - 0.34 
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erate (see TableO. This range corresponds to wavelengths from 
3.57 mm to 0.80 mm. The correlator can subdivide the frequency 
bands into a large number of spectral channels (~ 16 000) al- 
though lower spectral resolution can easily be achieved by re- 
configuring the correlator or by post-correlation spectral aver- 
aging of the data. A minimum integration time of 16 ms might 
be possible whereas switching between different bands will take 
~ 1.5 s. 

The field of view is defined by the primary beam size of an 
ALMA antenna, being 21" at A - 1 mm. That is sufficient to 
observe the interior of an internetwork region. In interferomet- 
ric mode ALMA will provide an angular resolution of 0'.'015 
to 1 '.'4, depending on antenna configuration. The angular reso- 
lution corresponds to ~ 10 km to ~ 1000 km on the Sun. The 
number of baselines is 1225 for the main array with 50 large 
antennae (see Eq. ([T)). With Eq. ^ we estimate the effective 
spatial resolution for images of maximum reliability, (Aor)^, to 
be of the order of ~ 0'.'42 at a wavelength of 1 mm. The resolu- 
tion for 0.3 mm and 3 mm would be 0'.'13 and 1'.'27, resp. (see 
also Table [TJ. A number of 64 antennae, as originally planned, 
would result in 2016 baselines. The total number of baselines, 
when including AC A, could be up to 2016 H- 120, reaching an 
effective resolution of ~ 0'.'32 at /I = 1 mm. 

2.2. CARMA 

The two millimetre arrays of Owens Valley Radio Observatory 
(OVRO) and of the Berkeley-Illinois-Maryland Association 
(BIMA) were merged to form the Combined Array for Research 
in Millim eter-wave Astronomy (CARMA) at Ce dar Flat, 
CaUfornia (I Woodv et al.l2004l:lBeaslev & Vogell2003 h. This het- 
erogeneous array consists of six 10.4 m antennae (OVRO) and 
ten 6.1m antennae (BIMA), resulting in 105 baselines. There are 
receivers available for 115 Ghz (2.6 mm) and 230 Ghz (1.3 mm). 
CARMA will be supplemented with a subarray consisting of 
8 3.5 m antennae from the Sunyaev-Zeldovich Array (SZA) 
for observations at 115 GHz but also 35 GHz (~ 1 cm) with a 
total of 276 possible baselines. Different array configurations 
with spacing from down to 5 m to 1.9 km are possible. An an- 
gular resolut ion of 0"! can be reached (for the 230 GH z A- 
array). See Wood v et all ( |2004|) . iBeaslev & Vogel ( l2003h . and 
http : //www . mmarray . org for more information. 

2.3. EVLA 

During the first project phase the Expanded Very Large Array 
(EVLA, see http://www.aoc.nrao.edu/evla) consists of 
the 27 VLA antennae with primary reflector diameters of 25 m, 
providing 351 baselines. In the second phase the array will be 
supplemented with eight new antennae, resulting in baselines of 
up to 350 km. The array will then have a very high angular reso- 
lution of 0'.'004 at /I = 6 mm. The accessible frequencies range 
from 1.0 GHz to 50 GHz, corresponding to 30 cm to 6 mm. The 
correlator will provide some thousand frequency channels. 

2.4. FASR 

The Frequency-Agile Solar Radiotelescope (FASR, see 
http : //www. ovsa.njit . edu/fasr) will combine three types 
of antenna, necessary to cover the large range from 30 MHz 
(A * 10 m) to 30 GHz (A ^ I cm). The sub-array for the high 
frequencies will consist of ~100 antennae of 2 m diameter 
each, setting up roughly 5000 baselines. The maximum antenna 



spacing will be 6 km. At 30 GHz an angular resolution of 
~0'.'66 can be reached. FASR will produce image sequences 
including polarisation information with a time resolution of 
less than 0.1 s. The construction is expected to be completed in 
2010. 



2.5. RAINBOW 

The six transportable lOm-antennae of the Nobeyama 
Millimeter Array have been linked to the local 
45 m-antenna to form the RAINBOW interferometer 
(http://www.nro.nao.ac.jp). The maximum baseline 
length is of the order of 400 m. RAINBOW can access 
wavelengths from 1.3 mm to 3.5 mm. 



3. Numerical models 

The numerical 3-D models used here are ca lculated with 
the radiation hydrodynamics code CO^BOLD jFrevtag et al.l 
200| i). Most of the s t udy re fers to the non-magnetic model 
by iWedemever et al.l ( |2004 hereafter W04, model A), 
whereas only o ne snapsho t is ta ken each from the mod- 
els by [Sc haffenb erger et alj (l2006l hereafter S06, model B) 
and iLeenaarts & Wedemever-Bohml (12006b ', hereafter LW06, 
model C). The models are described below (see also Table|2]i. 

For all models a grey, i.e., frequency-independent radiative 
transfer is used. The lateral boundary conditions are periodic in 
all variables, whereas the lower boundary is "open" in the sense 
that the fluid can freely flow in and out of the computational do- 
main under the condition of vanishing total mass flux. The spe- 
cific entropy of the inflowing mass is fixed to a value previously 
determined so as to yield solar radiative flux at the upper bound- 
ary. The upper boundary is "transmitting" in the hydrodynamics 
case (models A, C) but "closed" for the MHD model B, i.e., 
reflecting boundaries are applied to the vertical velocity, while 
stress-free conditions are in effect for the horizontal velocities. 

Please note that we define the term photosphere as the layer 
between km and 500 km in model coordinates, and the term 
chromosphere as the layer above. The origin of the geometric 
height scale is chosen to match the temporally and horizontally 
averaged Rosseland optical depth unity for each model individu- 
ally. The computational time step is around 0. 1 to 0.2 s in hydro- 
dynamics case (models A, C) and an order of magnitude smaller 
for the MHD model B. 

The utilised models still have short-comings concerning the 
energy balance of the chromosphere but nevertheless can give a 
first idea of the small-scale structure and dynamics of this layer. 
Refer to Sect. |6.1| for a discussion of the hmitations of the mod- 
elling. 



Table 2. Numerical models used for this study: publication de- 
scribing the individual models in more detail, short description, 
code used for intensity synthesis, number of time steps n,, and 
heliocentric angles (ju = cos 6). 



model 


pub. 


description 


synthesis 
code 


n, 




A 


W04 


non-magnetic 


LinforSD 


60 


0.2- 1.0 


B 


S06 


magnetic, 


LinforSD 


1 


1.0 






|Sol= lOG 








C 


LW06 


non-equilibrium 


RH 


1 


1.0 






H ionisation 
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Fig. 1. Horizontal maps from the non-magnetic model A (left column) and the somewhat smaller MHD model B (right column). 
Panels a) and f) show the gas temperature at a geometrical height of z = 1000 km. In panels b-e) the brightness temperature for 
model A is displayed for wavelengths of 0.3 mm, 1 mm, 3 mm, and 9 mm, respectively. The brightness temperature is also shown 
for model B at wavelengths of g) 0.3 mm, h) 1 mm, and i) 3 mm. The magnetic field strength \B\ at a height of z = 1000 km in 
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Fig. 2. Horizontal maps from the non-magnetic model C with electron densities treated in LTE (left column) and with the time- 
dependent non-equilibrium electron densities (NLTE, right column), which result from the simulation. The wavelengths from top 
to bottom are 0.3 mm, 1 mm, 3 mm, and 9 mm, respectively. Please note that the average intensity contribution function for panel d 
partly exceeds the upper boundary of the model. 



3.1. Field-free hydrodynamic model 

Model A (W04) does not include magnetic fields. The computa- 
tional domain consists of 140 x 140 grid cells in horizontal (x, y) 
and 200 in vertical direction. While the horizontal resolution is 
constantly 40 km, the vertical resolution varies from 46 km at the 
bottom in the upper convection zone at z = -1400 km to 12 km 
for all layers above z - -270 km. The top of the model is located 



at a height of z = 1710km in the middle chromosphere. The 
horizontal extent is 5600 km, corresponding to an angle of 7'.'7 
in ground-based observations. For the analysis presented here, 
we use a partial sequence with a time increment of 10 s and a 
duration of 600 s. The model chromosphere is characterised by 
a mesh-like pattern of hot shock fronts and cool post-shock re- 
gions inbetween (see Fig.[T^). 
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3.2. Magnetohydrodynamic model 

Model B (S06) includes a weak magnetic field with an aver- 
age flux density of 10 G. The computational domain is some- 
what smaller than for model A. It extends over a height range of 
2800 km of which 1400 km reach below the mean surface of op- 
tical depth unity (the same as for model A) and 1400 km above 
it. The horizontal extent is 4800 km x 4800 km. With 120^ grid 
cells, the spatial resolution in the horizontal direction is con- 
stantly 40 km, whereas it varies between 50 and 20 km in the 
vertical direction, thus accounting for the varying scale heights 
in the atmosphere. The initial model has a homogeneous, verti- 
cal, unipolar magnetic field with a flux density of lOG which is 
superposed on a previously computed, relaxed model of thermal 
convection very similar to model A. In the course of the sim- 
ulation the convective flows advect the magnetic field towards 
the intergranular lanes where stronger flux concentrations ("flux 
tubes") build up. In contrast, the model chromosphere is charac- 
terised by a more homogeneous but more rapidly varying field 
distribution (see Fig. [1]). The gas temperature still is similar to 
the one in model A (see Fig. [if). 

3.3. Model with non-equilibrium hydrogen ionisation 



Model C (LW06, see also iLeenaarts & Wedemever-Bohml 
l2006ah is identical to model A with respect to the numerical 
grid, the radiative transfer, and the hydrodynamics solver. It also 
does not contain magnetic fields. The only difference compared 
to model A is that the hydrogen ionisation is not treated in LTE 
but in non-equilibrium (non-LTE) by solving the time-dependent 
rate equations for a six level model atom with fixed radiative 
rates. Hydrogen is treated as minor species so far, i.e. the non- 
equilibrium hydrogen ionisation has no back-coupling on the 
equation of state and on the opacities. The gas temperature distri- 
bution is thus (statistically) the same as for model A. In contrast 
to model A, however, the code outputs the electron densities and 
hydrogen level populations for the LTE case and the NLTE case 
for each grid cell for model C. The electron density contribu- 
tion of hydrogen directly results from the non-LTE computation, 
whereas the contributions of other chemical species are treated 
in LTE (see LW06 for details). The electron density in the model 
chromosphere varies much less in non-LTE compared to LTE. 

4. Intensity synthesis 

Continuum radiation at (sub-)mm wavelengths is mainly due 
to thermal free-free emission and originates from the chromo- 
sphere and upper photosphere. The opacity is mostly due to 
free-free processes (interaction of ions and free electrons), in- 
cluding free-free (interaction of neutral hydrogen atoms and 
free electrons). Owing to the large wavelength and thus small 
frequency the condition 



(3) 



is fulfilled so t hat the Rayleigh- Jeans appro ximation can be used. 
According to 'Rv bicki & LightmanI (|2004 see Eq 5.19a) and 
[Mihalas (cf. 1978, p. 102), the opacity coefficient for thermal 
free-free bremsstrahlung (ion-electron interaction) can then be 
written as 



«e »i 
" ^ gas 



(4) 



ticular the dependence of Kg on and v"^ is of importance for 
the optical depth at wavelengths around 1 mm and thus of par- 
ticular interest for the study presented here. The free-free pro- 
cesses are due to collisions with electrons and depend on the lo- 
cal thermodynamic state of the electrons. The ratio of emission 
to absorption processes is thus in local thermodynamic equilib- 
rium (LTE), and the source function is Planckian. Another con- 
sequence of Eq.[3]is that the contribution to the emergent inten- 
sity / brightness temperature at (sub-)mm wavelengths is linearly 
related to the local gas temperature in the contributing height 
range. The electron densities, however, can deviate from the LTE 
values. That is caused by (i) ionisation by a non-Planckian radi- 
ation field, in particular in the Balmer continuum, and (ii) the 
long recombination timescales that hinder the hydrogen ionisa- 
tion degree to foll ow the faster d ynamic changes of the atmo- 
spheric conditions jCarlsson & St ein 2002, LW06). This devia- 
tion from equilibrium has an effect on the optical depth and thus 
on the effective formation height of the radiation via the absorp- 
tion coefficient (see Eq. (|4|). For this first quaUtative study, we 
neglect the resulting effect on the opacity for the intensity syn- 
thesis from model A and B but investigate the effect for model C 
(see Sects. |5^and|6l). 

For the radiative transfer calculations for model A 
and B we use LinforSD, a 3D LTE spectrum synthesis 
code developed by M. Steffen and H.-G. Ludwig, which 
is originally based on the Kiel code LINFOR/LINLTE 
(see http : //www . aip . de/~mst/lin£or SDjnain . html). 
Electron densities are calculated under the assumpt i on of LTE. 
For model C we use the RH code by lUitenbroekl (l2000h . We 
calculate a snapshot with LTE electron densities but also with 
the non-equilibrium electron densities, which are a direct result 
of the simulation with non-equilibrium hydrogen ionisation (see 
Sect.|l3]l. 

Continuum intensity images are calculated at the four wave- 
lengths 0.3mm (~ lOOOGHz), 1mm (~ 300GHz), 3mm (~ 
100 GHz), and 9 mm (~ 33 GHz). The heliocentric angle is here- 
after referred to as the inclination angle 9 and its corresponding 
cosine // - cos 6. For model A we consider five different posi- 
tions on the solar disk from yU = 1 .0 (disk-centre) to 0.2 (near 
limb) with an increment of AyU = 0.2. With a number of 60 snap- 
shots this results in a total of 60 x 4 x 5 intensity maps. For 
model B and also the LTE and NLTE case for model C only disk- 
centre maps for one snapshots are computed, giving 4 intensity 
maps in each case. 

For convenience, intensities I,\ in units of 
ergcm"^ s"' A"' sr"' (as output by LINFORSD) are con- 
verted to brightness temperature Tb via a relation derived from 
the Kirchhoff'-Planck function in the limit of Eq. ([3]l: 



A" 



■h- 



(5) 



Here, A, k-^, and c are the wavelength, the Boltzmann constant 
and the speed of light, respectively. The intensity contribution 
function C along a vertical view angle in - \ .0) is defined as 

/f(z)p(z)S,(z)e-^^«. 



C(z) 



(6) 



where and n\ are the number densities of electrons and ions, 
resp., and v and Tgas are frequency and gas temperature. In par- 



The above quantities are opacity k (extinction coefficient per 
mass unit), density p, Planck function S,;, and optical depth t,i 
which are functions of geometrical height z. Integration of the 
contribution functions over all heights yields the emergent in- 
tensity. 

Examples for the resulting intensity maps are shown in Fig.[T| 
for models A and B at disk-centre, in Fig.|2]for model C, and for 
different inclination angles for model A in Fig.|9l 
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d: Tu, A, = 3.0 mm 




e: Th, X, = 9.0 mm 




a:T„3s,z = 500 km, 1000 km 
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Fig. 3. Temperature distribution in model A: Histograms for a) 
gas temperature in horizontal slices at geometrical heights of 
z - 500 km (dotted) and z - 1000 km (solid), and b-e) bright- 
ness temperature at wavelengths of 0.3 mm, 1 mm, 3 mm, and 
9 mm, respectively (all disk-centre). All time steps of the anal- 
ysed 3-D model sequence are taken into account. The vertical 
dashed lines mark the positions of the individual hot and cool 
peaks. 



5. Results 

5.1. Brightness temperature distribution 

The intensity or brightness temperature maps for the different 
wavelengths all exhibit the complex pattern of hot/bright fila- 
mentary structures and cool/dark regions inbetween that is al- 
ready seen in the gas temperature cuts through the model chro- 
mospheres. See Fig. [T] for an example of a time step from 
model A and model B, and Fig. |2] for the LTE and NLTE cases 
from model C, respectively. Nevertheless, there are differences 
in the brightness temperature distribution which we quantify 
here by means of the horizontal and, in case of model A, tempo- 
ral average {T\,} and the relative intensity contrast 



01) 



(7) 



where ^rb.rms is the brightness temperature fluctuation. The val- 
ues are listed in Table [3] for all models, wavelengths, and incli- 
nation angles. As we shall see later in more detail, the longest 
wavelength is only of limited meaning in the present study since 
the formation height range partly exceeds the vertical extent of 
the numerical models. 

The average brightness temperature for model A increases 
with wavelength A which is due to the corresponding change 
in effective formation height (see Sect. 15.21 ). The same effect is 
found when going from disk-centre towards the limb, i.e. with 



increasing inclination angle 0. The LTE case for model C be- 
haves in a similar way whereas the average temperatures de- 
crease with wavelength for the corresponding NLTE case and 
also model B. 

The absolute brightness temperature fluctuation 5rb,i-ms tends 
to increase with wavelength at or close to disk-centre while this 
trend reverses closer to the limb. For model A this means an 
increase of the rms temperature amplitude from ~ 800 K to 
~ 11 00 - 1200K. The effect is more pronounced for models B 
and C. The relative brightness temperature contrast 6Th,ims/{Tb} 
depends in a similar way on wavelength and inclination angle. 
Despite the large temperature fluctuations in the layers appar- 
ently sampled by the long wavelengths (see Sect. I5.2l l. the con- 
trast reduces to only ~ 12 % for A = 9 mm close to the limb. 
In contrast, the maps from model A at disk-centre show con- 
trasts of up to 24 % with exception of the shortest wavelength 
that provides 19 % at maximum. Model B has much larger fluc- 
tuations with contrasts from 27 % to 34 % implying a detectable 
influence of the weak magnetic field on the gas and brightness 
temperature distribution. The contrast in the snapshot with LTE 
electron densities from model C is similar to model A, although 
with a larger spread with wavelength. This can be explained by 
the fact that the contrast for model A relies on 60 snapshots but 
only on a single one for model C. The trend of the contrast to 
increase with wavelengt is even more pronounced in the NLTE 
case. The more homogeneous electron densities reduce the con- 
trast at /I = 0.3 mm to only ~ 13 % but increase it for the other 
wavelengths. The different behaviour of the shortest wavelength 
is related to the fact that the bulk of emission originates pre- 
dominantly from layers where the amplitudes of the propagating 
shock waves are still relatively small (upper photosphere) - in 
contrast to the other wavelengths which sample higher layers 
(see Sect. l5.2l i. 

The less corrugated surface of equal optical depth results in 
less contributions from layers with higher gas temperature fluc- 
tuations. As a consequence the NLTE case provides a cleaner 
measure of the layer near the classical temperature minimum at 
the top of the model photosphere where the temperature fluc- 
tuations are small. At the longer wavelengths, only the shock- 
patterned chromosphere is sampled in the NLTE case and also 
in LTE there are significant contributions from that layer. In con- 
trast to the region of low amplitude around the classical temper- 
ature minimum the chromosphere shows a similar pattern of hot 
shock fronts and cool post-shock regions for all heights above 
z ~ 700 km. In LTE, the optical depth roughly follows the tem- 
perature gradients so that the chromosphere is sampled on cor- 
rugated surface of equal optical depth. The corresponding sur- 
faces in NLTE vary much less in height and thus cut through the 
temperature fluctuations instead of following them, explaining 
the higher contrast in NLTE compared to LTE. The difference 
in intensity distribution is shown more detailed in form of his- 
tograms in Fig. |3]for the four wavelengths (all for disk-centre 
only) and for gas temperature at geometric heights of z = 500 km 
and 1000 km (uppermost panel) for model A. As described in 
more detail by W04, the distribution of gas temperature exhibits 
a cool and a hot component with an intermediate range. The 
two components are due to a cool background and hot shock 
fronts. W04 also provide histograms for other heights (see Fig. 7 
therein). The histograms for continuum intensity show two com- 
ponents, too, although with varying amplitudes. At A - 0.3 mm 
a strong low-intensity component is visible whereas it is hard to 
define a high-value peak at all. The distribution for A - I mm 
is most similar to the gas temperature histogram at z = 1000 km 
with peaks in almost the right proportion. In contrast, the high- 
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Table 3. Average brightness temperature (Tb), rms variation dT^ i-^s, and brightness temperature contrast ^rb,ims/(7'b) for all wave- 
lengths A and disk-positions fi derived from all snapshots. Due to the linear conversion at given A the values can also be interpreted 
as intensity contrast dli-ms I {!)■ In some cases the corresponding average contribution function partly exceeds the upper boundary of 
the model. The missing contribution is in the range of 0.5 % to 2 % for data marked with * and > 2 % for data marked with **. 
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A [mm] 
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Fig. 4. Normalised contribution functions for continuum inten- 
sity on the geometrical height scale calculated with Linfor3D 
from the model by W04 (average over all horizontal positions 
and time steps in the analysed data sample) at wavelengths of 
a) 0.3 mm, b) 1 mm, c) 3 mm, and d) 9 mm for ju = 1 .0 (thick 
soHd), |/ - 0.8 (dotted), yu = 0.6 (dashed), yi - 0.4 (dot-dashed), 
yu = 0.2 (triple-dot-dashed). The assumption of LTE was made 
for the calculations. 



intensity peak is more pronounced than the low-value compo- 
nent in case of the longer wavelengths (/I = 3 mm, 9 mm). The 
differences between the wavelengths can be understood if one 
considers the different height ranges that contribute to the con- 
tinuum intensity (see Sect. |5.2| i. The low-value peak in gas tem- 
perature at z = 1000 km is at lower values than for the bright- 
ness temperature distributions. A linear conversion between both 




200 400 600 800 1000 1200 1400 1600 
z[km] 



Fig. 5. Normalised contribution functions for continuum inten- 
sity on the geometrical height scale for average from model A 
(see Fig. ID here dotted line), for the snapshot from MHD 
model B (dashed), and for the snapshot from the model with 
non-equilibrium hydrogen ionisation. For the latter the electron 
densities were first calculated with RH under the assumption of 
LTE (dot-dashed) and then with non-equilibrium electron densi- 
ties available from the model (thick solid). The panels show dif- 
ferent wavelengths: a) 0.3 mm, b) 1 mm, c) 3 mm, and d) 9 mm. 



temperatures due to the validity of LTE (see Sect.|4]l should pro- 
duce identical distributions for gas temperature and coiTespond- 
ing brightness temperature at the same height. The differences in 
Fig. |3] however, illustrate the influence of the extended forma- 
tion height ranges (see Sect. l5.2l i which effectively mixes contri- 
butions from different layers instead of sampling the temperature 
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Fig. 6. Formation heights Zcc as function of brightness temper- 
ature for the snapshot from model C with LTE electron densi- 
ties (left column) and non-equilibrium electron densities (right 
column). The range in brightness temperature and formation 
height covered by all horizontal positions was divided into bins 
of Az - 12 km and AT^ = 50 K. The number of spatial positions 
with values within the same bin (i.e., the density function) is 
shown as grey-scale for the wavelengths of 0.3 mm (a,e), 1 mm 
(b,f), 3 mm (c,g), and 9 mm (d,h). The solid lines represent the 
average brightness temperatures and formation heights, whereas 
the dotted hnes mark the 1 cr deviation. 



at a fixed geometrical height, as it is the case for the displayed 
gas temperature maps. 

5.2. Contribution functions and formation heights 

In the following we describe the efifective formation heights by 
means of Zmaxc, the height of the maximum of the horizontally 
(and temporally) averaged contribution functions, and by Zcc, 
the height of the centroid of the (spatially resolved) contribution 
functions, where Zcc is defined as 



/z'-C(z')*' 
/ C{z')dz' 



(8) 



All available horizontal positions (and time steps) are used for 
the calculation of the average (zcc) and the variation 5(Zcc)ims- 
The results are summarised in TableH) The horizontally and tem- 
porally averaged contribution functions for model A are shown 
in Fig. |4] for all wavelengths and inclination angels. For disk- 
centre (ji - 1.0), the continuum intensity at 0.3 mm originates 



mainly from layers near the classical temperature minimum re- 
gion with a contribution peak at a height of Zmaxc = 48 1 km. 
Additionally, there are small contributions from the low chromo- 
sphere. The picture is in principle the same for a wavelength of 
1 mm, although the peak is located somewhat higher at a height 
of 588 km and the chromospheric contribution is larger. Hence, 
these both wavelengths map mainly the thermal structure of the 
top of the photosphere and the low chromosphere. Next to a 
maximum at Zmaxc = 695 km the emission at /} = 3 mm origi- 
nates also from an extended height range in the chromosphere 
with a subtle secondary maximum at 1120 km. At A - 9 mm 
all chromospheric heights contribute almost equally to the inten- 
sity with a maximum at 1348 km and a less pronounced, slightly 
smaller peak (0.98) at a height of 789 km which might be com- 
pared to the peaks for the shorter wavelengths. The contribu- 
tion function implies that there would still be significant emis- 
sion from layers above the computational domain. Therefore 
the intensity images for A - 9 mm, as shown in Fig. [T] and 
also the results derived for this wavelength are incomplete to 
some extent and thus of only limited significance. According to 
ILoukitcheva et al. (2004), emission at /I > 8 mm does even origi- 
nate from the transition region which is not included in the mod- 
els used here. Consequently, a more extended model would be 
necessary for correctly modelling emission at long wavelengths. 

Owing to the complicated thermal structure and the corre- 
sponding LTE electron density, a large height range is involved 
in emitting radiation at a certain wavelength, causing the contri- 
bution function to be rather complicated. The peaks of the indi- 
vidual low-altitude components of those functions for model A 
(see Fig.lHi move upwards by roughly 100km for each increase 
of factor 3 in wavelength. But additionally the high-altitude 
contribution, say above ~ 900 km, grows strongly with wave- 
length. The mean formation height (zcc) therefore increases by 
~ 200 km with each factor 3 in wavelength. 

The average contribution functions for an inclined view 
in < 1.0) in model A are shown in Fig. 21 too. At all wave- 
lengths a similar behaviour is obvious: The low-altitude peak (as 
most clearly seen for short wavelengths at large /i) moves higher 
up and at the same time decreases its relative contribution un- 
til the high-altitude range prevails. Although the latter is much 
broader, in contrast to the relatively sharp low-altitude peak, a 
maximum can be defined in most cases. The height of this max- 
imum increases with decreasing yU, just as for the low-altitude 
contribution. The mean formation height (zcc) increases by ap- 
proximately 200 km with each factor 3 in wavelength for all in- 
clination angles. 

The contribution functions for the snapshot from model B, 
which are shown in Fig.|5] are similar to those for model A, con- 
sidering that one compares an average over 60 time steps with a 
single snapshot. Also the formation heights agree generally (see 
Table |4|i, although a factor 3 in wavelength only increases the 
mean formation height by A(zcc) ~ 130 km (Az,naxc ~ 140km). 

Also the contribution functions for the snapshot from 
model C are similar when using LTE electron densities for the 
intensity synthesis (see dot-dashed line in Fig. |5]l. Using the 
non-equilibrium electron densities from the simulation, how- 
ever, produces different contribution functions for wavelengths 
of 1 mm and longer. At A - 0.3 mm LTE and non-LTE both are 
very similar to the other models, reflecting the fact that LTE is 
still a reasonable assumption for electron densities at that height 
(LW06). At the longer wavelengths the relevant height ranges 
of the contributions functions seem to be more localized in the 
non-LTE case and do not exhibit such a broad range and double- 
peaked shape as for the LTE case. The difference is caused by 
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Table 4. Heights of maximum contribution z^axc (see Fig. |4]i, average and variation of the heights of the centroids of the spatially 
resolved contribution functions, (zcc) and (5(zcc)rms, resp., for all wavelengths A and disk-positions fi derived from all snapshots. The 
height values in square brackets only represent secondary maxima. In some cases the corresponding average contribution function 
partly exceeds the upper boundary of the model. The missing contribution is in the range of 0.5 % to 2 % for data marked with * and 
> 2 % for data marked with ** . 
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the electron densities which have direct influence on the optical 
depth. While in LTE the hydrogen ionisation and with it the elec- 
tron density follow the strong temperature variations of the prop- 
agating shock waves, they do vary much less in non-LTE and 
tend to be at values set by the conditions in the shocks (LW06). 
The effect is that surfaces of constant optical depth at wave- 
lengths around 1 mm are much less corrugated. Consequently, 
a much smaller height range is sampled in non-LTE whereas in 
LTE the strongly varying optical depth mixes contributions from 
low and high layers, depending on the thermal structure along 
the hn e of sight (see Fig. 3 by Leenaarts & Wedemever-Bohml 

The formation heights Zccix,y) of all columns in model C 
are shown in Fig. |6]in dependence of the resulting brightness 
temperature Ti,{x,y). The strong variation of the electron density 
and thus the optical depth in the LTE case makes the formation 
height range increase with brightness temperature. This trend is 
clearly visible for all wavelengths in the left column of Fig. |6] 
(LTE). A high brightness temperature is connected to high gas 
temperature (in a shock front) and thus to an increased opacity. 
Essentially the atmosphere gets optically thick already high up 
in the atmosphere if a shock wave is in the line of sight, whereas 
the cooler temperatures in the post-shock regions allow to look 
at much deeper layers. In NLTE (right column) this effect is less 
pronounced for the short wavelengths and essentially absent for 
the long wavelengths as the fluctuations of the non-equilibrium 
electron density are smaller The spread in height is reduced in 
the NLTE case compared to LTE, e.g. from i5(zcc)rms = 197 km to 
only 1 17 km at a wavelength of 1 mm. Still the variation prevents 
a strict correlation between wavelength and corresponding for- 
mation height that would be highly desirable for the derivation 
of the height-dependent three-dimensional thermal structure of 
the atmosphere from observations. Nevertheless, the (more real- 
istic) NLTE calculations for model C imply it could be possible 
at least in a statistical sense. The wavelength dependence of the 
average formation height remains very similar to the results from 
model A. An increase by a factor of 3 in wavelength results in 
an increase of A{zcc) ~ 190 km (Az^axc ~ 220 km). 

The clear tendency of sampling higher layers with increas- 
ing wavelength is expected as the opacity goes with the wave- 
length squared. In principle, the sampled height can be related to 
the density stratification. An increase of factor 3 in wavelength 



would then shift optical depth unity by that height after which 
the density decreased by 3^. In the models this height is mostly 
between 220 km (low chromosphere) and 400 km (towards high 
chromosphere) and thus agrees well with the values directly de- 
rived from the average contribution functions. 

5.3. Temporal evolution 

The analysis of the temporal behaviour relies on model A as 
only single snapshots are available for the other models. The 
dynamics of model C are the same and also the MHD model B 
is very similar in this respect. 

Pattern evolution time scale The temporal evolution of the in- 
tensity pattern can be quantified in terms of a time in which the 
autocorrelation decays to the fraction 1/e as it has been done in 
W04 for gas temperature. Here, we apply the same procedure 
to time sequences of brightness temperature for all wavelengths 
and all inclination angles from model A. The time scales are 
mostly around 23 s to 24 s for /u > 0.2 and go down to 19 s for 
fi = 0.2 and A = 3 mm. The extreme case of 17.4 s for ju - 0.2 
and A - 9 mm must regarded as uncertain as it is influenced 
by possible artefacts due the only partially included formation 
height range. The gas temperature in the model chromosphere it- 
self exhibits time scales of the same order (20 s to 25 s) and stays 
almost constant throughout the chromosphere with a tendency 
towards higher values for the lower layers (see W04). Hence, it 
is not surprising that the intensity image sequences reproduce 
roughly the same time scales for the different wavelengths de- 
spite the different contribution functions. 

Brightness temperature variations The dynamic nature of the 
model chromosphere is very obvious if one looks at the tem- 
poral variation of the brightness temperature in Fig. [T] (see also 
Fig. 10 of W04 for the chromospheric gas temperature varia- 
tion). The amplitude of the fluctuations, expressed quantitatively 
by means of Tb^rms (see Table [3] compare with Fig. 9 in W04), 
is generally smaller in the lower layers than further up, since the 
shock waves steepen on their way upwards into the thinner atmo- 
spheric layers. The propagating wave fronts are clearly visible as 
bright streaks in the lower panel which shows the gas tempera- 
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Fig. 9. Variation of continuum intensity maps 
at a wavelength of 1 mm with inclination angle 
(fi = cos d). The thick circles with solid dot in 
the sketches above each panel indicate the po- 
sitions on the solar disk. The different panels 
sizes are due to foreshortening. 
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Fig. 7. Temporal variation of excess brightness temperature at 
a chosen horizontal position (x - 4460 km, y - 2020 km) in 
model A for all four wavelengths (a-d) at disk-centre. In each 
panel the standard deviation is noted. The lower panel (e) shows 
the gas temperature (colour legend below) in the selected col- 
umn as function of time and height. The vertical long-dashed 
lines mark shock fronts at A - 0.3 mm and (together with the 
short-dashed lines at equi-distant times) help to identify the up- 
ward propagating fronts after a short delay in the brightness tem- 
perature for the longer wavelengths. 
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Fig. 8. Gas temperature profiles along the horizontal x axis (up- 
per panel) and the vertical z axis for the same snapshot as in 
Fig. [T]at (x,y,z) = (4460km, 2020km, 1000km) (solid) and . 
(x,y,z) = (5220km, 2020km, 1000km) (dashed). The positions 
are marked with dotted lines. The arrows at the right indicate the 
line of sight direction at the limb (fi = 0.0) and at disk-centre 
(ju = 1-0). 



The vertical dashed lines in the figure mark shock fronts at 
A = 0.3 mm, which are formed low in the atmosphere (see lower 
panel). The longer wavelengths show an excess in brightness 
temperature shortly afterwards. The time difference with respect 
to 0.3 mm increases with wavelength and thus with formation 
height, clearly indicating upward propagating shock waves. 

The vertical gas temperature stratification for the selected 
column is shown in Fig. [8] at a time of 480 s. The shock wave 
at a height of z = 1000 km appears as brightening at a wave- 
length of 1 mm at f = 480 s (Fig. |2ll and shortly afterwards at 
the longer wavelengths. The absolute temperature amplitudes of 
the shock waves, however, depend on details of the numerical 
modelling of the chromosphere (see discussion in Sect. 16. 11 1. The 
determination of brightness temperatures from observations will 
thus provide an important test for the numerical models. 



ture in the selected column as a function of height and time. The 
other temperature enhancements are caused by interaction with 
neighbouring wave fronts and the fact that the waves do not only 
move vertically (inside the column) but also laterally and thus 
can leave and enter the selected column sideways. 

Since the continuum intensity at A = 0.3 mm emerges from 
the high photosphere and low chromosphere, it exhibits smaller 
variations of the related brightness temperature compared to the 
remaining wavelengths which sample higher layers. For the par- 
ticular example in Fig.|7]we find a standard deviation of ~ 640 K 
at A = 0.3 mm and around 1000 K for the longer wavelengths. 



5.4. Centre-to-limb variation 

In Fig. |9]a sequence of intensity imagesfor A - I mm with in- 
creasing inclination angle is shown, illustrating the changes as 
one goes from disk-centre towards close to the limb. At disk- 
centre (fi - 1.0) the chromospheric small-scale structure with 
its bright filaments and dark intermediate regions is best visi- 
ble, resulting in the highest intensity contrast (see Table [3]). But 
towards the limb the contrast decreases since the dark regions 
become less dominant. This can be explained if one considers 
that at disk-centre the chromosphere is seen from above and at 
(or close to) the limb from the side, making visible different ge- 
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ometrical aspects of the thermal structure of the model chromo- 
sphere. From above the model chromosphere appears as alter- 
nating pattern of hot and cool gas (see Fig. |8]i and thus bright 
and dark in continuum intensity. There are usually only one or 
two (or none) shock fronts in the line-of-sight, allowing to look 
deep into the intermediate dark regions. In contrast from the side 
(small jj) many "propagation channels" of shock waves overlap 
along the line-of-sight, resulting in much lower intensity con- 
trast (see Fig.|8]l. Observations of the centre-to-limb variation of 
the intensity pattern may thus give valuable information on the 
three-dimensional structure of the solar chromosphere. 




3 : -I D 



Fig, 10. Intensity image aX. A - \ mm from a single snapshot of 
model A at different spatial resolutions: a: 0'.'06 (size of com- 
putational grid cells in the original model), b: (X'S, c: 0'.'6, d: 
(y!9. The value range is the same for all panels. The brightness 
temperature contrast is noted right next to each panel. 



5.5. Spatial resolution 

A major disadvantage of past and present instruments for the 
(sub-) millimetre range is an only poor spatial resolution which 
hindered a detailed study of the small-scale structure of the solar 
atmosphere. In order to estimate the spatial resolution which is 
needed to resolve the pattern clearly present in the model chro- 



mospheres used here, we choose a representative intensity image 
from model A for A - I mm and degrade it artificially by convo- 
lution with a Gaussian kernel, mimicking instrumental proper- 
ties. The result is shown in Fig.[TO]for different widths (FWHM) 
of the Gaussian. The finest structures are akeady hardly visi- 
ble at a resolution of 0'.'3 whereas the larger mesh-like pattern 
is still discernible at 0'.'9. The brightness temperature contrast, 
which is noted next to the panels in Fig.[TOl reduces with spatial 
resolution Aa. This trend can be fitted with 

where D is the characteristic length scale of the chromospheric 
pattern (derived from the fit). This dependence is similar for all 
wavelengths and inclination angles although the length scale for 
the best fit is different. Except for fi = 0.2, where D x 290 km for 
all wavelengths, the length scale increases with wavelength. At 
disk-centre it increases by ~ 100 km for a factor three in wave- 
length from 730 km at /I = 0.3 mm to 1030 km at A - 9 mm. This 
can be explained by the fact that the propagating shock waves, 
which produce the pattern, only become visible at the bottom of 
the chromosphere, whereas in the higher layers the interaction of 
neighbouring wave fronts developed into a more distinctive pat- 
tern of larger mesh-size (see Fig. 2 by W04). On the other hand, 
D decreases with increasing inclination angle as one looks at the 
mesh-like pattern from above at disk-centre but sideways on the 
shock fronts close to the limb (see Sect. 15. 41 and Fig. |9]i 

6. Discussion 

6. 1 . Numerical modelling 

The results of this study demonstrate the constraints on spa- 
tial and temporal resolution for observations of the solar inter- 
network chromosphere with future (sub-)millimetre instruments 
but should be regarded qualitatively only since some important 
ingredients are still missing in the models used here. The ex- 
ample of model B implies that even weak magnetic fields, as 
they are expected for solar internetwork regions, will have most 
likely important (indirect) influence on structure and dynamics 
of the chromosphere. Nevertheless, the qualitative picture of pro- 
nounced inhomogeneities clearly present in the non-magnetic 
model A are also found for the weakly magnetic model B (see 
Sect. 16^211. 

A more severe limitation concerns the assumption of LTE 
for the radiative transfer and the equation of state, which is made 
here in order to keep the simulations feasible. At chromospheric 
heights the assumption is no longer valid. Instead a by far more 
involved treatment of non-equilibrium and radiative non-LTE ef- 
fects is required (e.g. hydrogen ionisation). This will certainly 
influence the energy balance and with it the thermal structure 
of the chromosphere. Taking into account the time-dependence 
of hydrogen ionisation, which is not instantaneous at chromo- 
spheric heights but rather is governed by long recombination 
time scales, leads to deviations from ionisation equilibrium and 
thus electron densities that can differ significantly from the cor- 
responding equilibrium values (Carlsson & Stein 2002, LW06). 
The effect on the optical depth at millimetre wavelengths is sig- 
nificant. The electron density affects the optical depth and thus 
determines which layer is effectively mapped. From the analysis 
of a recent simulations that accounts for non-equilibrium hydro- 
gen ionisation (model C, see LW06) we learn that (i) the height 
of optical depth unity at a wavelength of 1 mm is on average 
close to the LTE case and does not vary that strongly, (ii) the 
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intensity fluctuations and the average brightness temperature are 
reduced, (iii) the resulting intensity maps are still qualitatively 
similar (see Sect. 15. 11 1. However, the back-coupling of the non- 
equilibrium on the equation of state and on the opacities - a final 
step that has to be done - will most likely increase the peak tem- 
peratures of the shock waves. 

6.2. Magnetic field 

The snapshot from magnetohydrodynamic model B gives a first 
idea about the influence of the weak chromospheric magnetic 
field when compared to the non-magnetic model A (for A < 
9 mm). Already the brightness temperature maps in Fig.[T]show 
a similar pattern for both models. The characteristic length scale 
of the mesh-like brightenings, which are due to hot shock fronts, 
seems to be larger in the MHD case. The formation heights are 
very similar and are only slightly higher in model B although 
this might be due to the particular snapshot chosen. The aver- 
age brightness temperature is higher in model B for A = 0.3 mm 
and 1 mm but is the same than for model A at 3 mm. The av- 
erage actually decreases with wavelength for model B whereas 
it increases for model A. The rms brightness temperature vari- 
ations increase with A in both cases but are 400- 500 K higher 
in the MHD model. The different behaviour of the two models 
might be explained with the finding that the magnetic field - al- 
though still quite weak - becomes more important with increas- 
ing height (S06) and might influence the temperature distribution 
by suppressing power of the propagating waves. Also the field 
concentrations in the photosphere might change the wave exci- 
tation which should also become visible in the higher layers. The 
comparison of magnetic and non-magnetic models should be re- 
peated in more detail and with a larger number of snapshots in 
the future, but akeadya this first qualitative study suggests that 
high-resolution (sub-)millimeter observations can give (indirect) 
information on the magnetic field topology of the solar chromo- 
sphere. 

6.3. Syntlietic observations 

The computation of synthetic images that can be directly com- 
pared to observations requires more than a simple convolution 
with a Gaussian (Sect. |53] l. which is strictly speaking only valid 
in the limit of full M-v-coverage (see Sect.|2]i. For more realistic 
maps one needs to convolve a synthetic image with the primary 
beam response of the antenna used, Fourier-transform the result, 
and multiply with the M-v-coverage of the array configuration. 
The Fourier back-transform yields then the "dirty map". Such a 
detailed image construction, of course, depends on specific in- 
strumental properties and array configuration at the time of the 
observation. Nevertheless, synthetic images can be easily gener- 
ated to match the properties of a particular observation. 

6.4. Comparison to observations and otiier modeis. 

iLoukitcheva et al.l ('2004') compare the dynamic ID simulation 
by CS and the semi-empirical FAL model A with a large num- 
ber of observations (see references therein). They conclude that 
the FAL model A could also account for the measured bright- 
ness temperatures, at least when combined with the CS model. 
Although the brightness temperatures calculated from the FAL A 
stratification are still within the error bars of the observations, the 
values tend to be at the upper limit if not even exceeding it. The 
CS model matches the observations much better as it predicts 



much lower temperatures than FAL A. The discrepancy is still 
small for short wavelengths (A ~ 0.1 mm), where the intensity 
originates from deep in the atmosphere, but increase to the order 
of 3000 K for the longer wavelengths. 

The average brightness temperatures in our models is similar 
to the values for the CS simulation and the observations for the 
shorter wavelengths {A < 1 mm). For the longer wavelengths, 
however, CS, FAL A, and also the observed temperatures in- 
crease strongly with wavelength whereas our models yield con- 
siderably smaller average values of the order 4000 K - 5000 K. 
In case of A = 9 mm the discrepancy can partly be attributed to 
the fact that the continuum forming layers exceed the vertical 
extent of model A and B and also C in the LTE case although 
this argument does not hold for the NLTE case for model C. 
The transition region and its steep temperature gradient, which 
is included in the CS and FAL A models but not in ours, is a 
more likely source of discrepancy. In the CS simulation the tran- 
sition region starts around a height of 1500 km from where in our 
models a significant part of the emergent continuum intensity at 
long wavelengths originates. Consequently the temperatures at 
long wavelengths would be too low and increase of the average 
with wavelength not as steep as found by CS and FAL. Last but 
not least the brightness temperatures distribution is affected by 
uncertainties in the chromospheric energy balance due to a too 
simple radiative transfer and the still missing back-coupling of 
the non-equilibrium hydrogen ionisation to the equation of state. 
The latter would lead to higher shock peak temperatures. 

The brightness temperature variations derived in this work 
(see Table O tend to increase with waveleng th (at disk-centre) 
and thus show the same tendency as found bv Louki tcheva et al] 
(2004). They also explain this behaviour with different forma- 
tion height ranges increasing with wavelength. For the very long 
wavelengths (/i > 8 mm), however, the temperature variations in 
the CS model decrease agai n. The rms variation, shown in Fig. 4 
of Loukitchev a et al.l (|2006') as function of wavelength, actually 
has a maximum at around 2 mm with (Jrb.rms = 930 K. The cor- 
responding variations in the models presented here are signifi- 
cantly larger for all considered wavelengths. This is related to 
the differences in the profile and propagation of the shock waves 
in the ID model CS and our 3D models. While sawtooth-like 
shock profiles are clearly present in th e CS model (see Fig. 3 by 
ILoukitcheva et al.Ll20o4 and Fig. 1 bv lLoukitcheva et al.L 120061) . 
they are not so easily visible in the models used here (see Fig.|7]l. 
In 3D - in contrast to a ID model - shock waves propagate in 
all spatial directions and interact with neighbouring wave fronts 
and this way obscure the shock signatures in vertical columns. 
Nevertheless, there are excess temperatures of more than lOOOK 
in both kinds of model. The case of 0.05 mm for the CS simula- 
tion varies less and resembles more a regular oscillatory behav- 
ior as it originates from further down in the photosphere. 

6.5. Formation iieigiits 

Double-peaked c ontribu tion functions were found by 
ILoukitcheva et al.l (|2004|) for simulated brightness temper- 
ature at millimetre wavelengths for single snapshots from the 
ID simulation by CS (see Fig. 5 therein). A secondary maximum 
is also visible in the average contribution for A - \ mm. This 
agrees well with the contribution functions for the models used 
here. Secondary maxima are more obvious in the contribution 
functions for individual columns, whereas the functions shown 
in Figs. ID and |5] are smoother as they represent averages from 
3D calculations with a large number of columns. More than 
one peak appears if there is more than one shock wave in the 



14 



Sven Wedemeyer-Bohm et al.: Inter-network regions of the Sun at millimetre wavelengths 



line of sight or at least a shock wave above, i.e. at smaller 
optical depth, the average formation height range set by the 
background stratification. Multiple contribution peaks are thus 
a direct imprint of the highly dynamic simulations and details 
certainly depend on the properties of the shock waves. 

A comparison of the average contribution functions for the 
models by CS and FAL A (Figs. 6 and 7 in Loukitcheva et al.l 
[2004, resp.) shows that the formation heights are in line with 
our results with exception of A = 10 mm in FAL A which is 
formed at a height of ~ 2250 km in the transition region. In 
contrast we find good agreement for the longest wavelength 
in our models with the CS model. For the other wavelengths 
(0.3 mm - 3 mm), however, the average formation height tends 
to be slightly lower in our model s com pared to CS and FAL A. 
According to iLoukitcheva et al.l (1200 4). the emission at A > 
8 mm does originate from the transition region which is not in- 
cluded in the models used here. Consequently, a more extended 
model would be necessary for coiTectly modelling emission at 
long wavelengths. 

The formation height range for model A seems to be too 
broad to derive information about the thermal structure from 
observations. The situation improves when taking into account 
the non-equilibrium electron density contribution due to hy- 
drogen as this flattens the surface of equal optical depth. As 
the electrons in chromosphere are predominantly due to hy- 
drogen, non-equilibrium modelling of other electron donators 
(which are cuiTently treated in LTE) are thus expected to have 
only a minor influence on the electron number density. But 
even with an advanced non-equilibrium modelling the contri- 
bution functions extend over a significant height range, caus- 
ing the eff'ective formation heights to vary over ~ 100 km 
(m odel C, see TablelH). This is illustrated very clearly in Fig. 3 
by lLeenaarts & Wedemeyer-Bohml(l2006bh . 

The resulting brightness temperatures should thus be inter- 
preted very carefully as they represent the integrated physical 
state of an extended height range. Hence, a brightness tempera- 
ture can deviate significantly from gas temperature at a geomet- 
ric height in the case of an inhomogeneous and rapidly evolving 
chromosphere. 

Nevertheless, our study is promising in the sense that simul- 
taneous observations at many millimeter wavelengths can still 
serve as a statistical tomography of the chromospheric layers. 
Not only the wavelength dependence but also the centre-to-limb 
variation of the formation height might prove useful in this re- 
spect. 

6.6. Suggested observations 

Based on the results concerning the spatial resolution and effec- 
tive formation heights discussed above, we can suggest which 
will be the most promising instruments (see Sect. |2]i for observ- 
ing the small-scale structure of the solar chromosphere. 

Although RAINBOW does access the wavelength range of 
interest here, the small number of antennae will not allow for the 
high spatial resolution required for imaging the chromospheric 
fine-structure. For the same reasons CARMA cannot provide 
brightness temperature maps that have high spatial and tempo- 
ral resolution at the same time. EVLA does consist of more an- 
tennae, resulting in a much larger number of baselines and thus 
a much higher spatial resolution. The minimim wavelength of 
7 mm, however, is at the limit for our application as a signifi- 
cant part of the continuum radiation at this wavelength and be- 
yond originates from layers that are not or only partly within 
the computational domain of our models. Nevertheless, EVLA 



is certainly of interest for observing the upper chromosphere / 
transition region. In this respect it connects to the wavelength 
range covered by ALMA. The shortest accessible wavelength of 
FASR will be even 1 cm. Just as EVLA, it won't be able to map 
the low and middle chromosphere but will be of great value for 
investigating the solar corona at even higher spatial resolution 
flian EVLA. 

The by far most promising tool for imaging the chromo- 
spheric fine-structure at high cadence is ALMA as it can ac- 
cess the spatial and temporal scales implied by our study. The 
effective spatial resolution, however, might be critical for the 
longer wavelengths. A possible solution could be the application 
of Multi-Frequency Synthesis (MFS, see Sect.|2]i. 

The wavelength range of the individual bands of ALMA 
(see Sec. 12. It coiTesponds to an expected difference in forma- 
tion height range of ~ 100 km between lower and upper limit of 
the band. The thousands of spectral channels by ALMA, which 
are output simultaneously for each wavelength / frequency band, 
thus correspond to different layers, stacked with slightly dif- 
fering formation heights. Tomographic reconstruction methods 
might then be applied to derive the three-dimensional structure 
of the layer sampled by the band. 

Simultaneous observations in different frequency bands will 
not be possible but switching between bands can be as fast as 
1.5 s which is fast enough in case of the solar chromosphere. As 
upwards propagating waves show up in the different bands sub- 
sequently with increasing wavelength (see Sect. 15. 3t , it would 
be best to cycle through the different bands from high to low 
frequencies. Tracking waves through the atmosphere should be 
possible with ALMA. The subsequent appearance of particular 
features at the different wavelengths will enable the determina- 
tion of the propagation speed and might - in combination with 
measured amplitudes - give further insight into wave dissipation 
and with that into the acoustic heating contribution. Combined 
simultaneous observation runs with optical diagnostics for the 
low and middle photosphere have the potential to reveal the ex- 
citation sources of waves. 



7. Conclusions 

Based on recent three-dimensional (magneto-)hydrodynamic 
simulations we conclude that future (sub-)millimetre instru- 
ments like ALMA will be capable of linearly mapping the gas 
temperature distribution in the chromosphere. The temporal and 
spatial resolution of ALMA will be sufficient to resolve structure 
and dynamics on scales inherent in the models presented here. 

The inhomogeneous and dynamic nature of the (model) 
chromosphere causes the intensity contribution functions and 
with that the formation height ranges to be extended and to vary 
in time and space. The variations, however, stay moderate as the 
electron density, which mainly defines the optical depth at wave- 
lengths around 1 mm, is connected to the rather well-stratified 
ionisation degree of hydrogen. The formation height range in- 
creases with increasing wavelength. The same effect is found for 
decreasing inclination angle, i.e. higher layers are sampled close 
to the solar limb than at disk-centre. This behaviour could be 
exploited for a (statistical) tomography of the thermal structure 
of the solar atmosphere. In the case of ALMA the upper photo- 
sphere to middle chromosphere will be mapped with the differ- 
ent wavelength bands with high spatial and temporal resolution. 
The large number of spectral channels within each wavelength 
band will result in many simultaneous maps which correspond 
to continuously stacked formation heights. Tomographic recon- 
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struction techniques might then allow volume imaging of the so- 
lar atmosphere. 

Hence, it should even be possible to track propagating hydro- 
dynamical waves in the solar chromosphere and to study chro- 
mospheric oscillations in more detail. This will have import im- 
plications for the discussion of the heating mechanism not only 
for the Sun but for stellar chromospheres in general. The range 
of application is broad and does include the study of the fine- 
structure of the magnetic network, too, direct measurement of 
the magnetic field topology, and even full-disk mosaics. Thus, 
ALMA will be of great value for finally resolving the discussion 
about the dynamics and the true structure of the solar chromo- 
sphere. 
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